%% Numerical Integration of ODE

%Here we are using the Euler Method
web('https://en.wikipedia.org/wiki/Euler_method')

%Setup our problem

%First order ODE
lambda=2;
Fp=@(F)-lambda*F; %this is the derivative of F(t) at F

%Integration time and time step
T=10;
Dt=.1;
T0=0;
Tvec=(T0:Dt:T) %time vector

%Intialize results array
F=zeros(length(Tvec),1);

%Set initial values
F0=10;
F(1)=F0;

%This is the integration 
for tt=2:length(Tvec)
    F(tt)=F(tt-1)+Fp(F(tt-1))*Dt;
end

%Show our results and the known exact result. Try changing Dt. 
figure
plot(Tvec,F)
xlabel('T')
ylabel('F(t)')

ExactT=linspace(T0,T,10000)
ExactF=10*exp(-lambda*ExactT)
hold on
plot(ExactT,ExactF)
legend('Numerical','Exact')

%% Second Order ODE

T=20;
Dt=.01;

%This is the spring equation
k=1;
m=1;
vp=@(x)-(k/m)*x;
xp=@(v) v;

%Setup
T0=0;
Tvec=(T0:Dt:T)
X=zeros(length(Tvec),1);
V=zeros(length(Tvec),1);

X0=10;
V0=0;

X(1)=X0;
V(1)=V0;

%This is the integration.  In this case, we update both velocity and
%position
for tt=2:length(Tvec)
    X(tt)=X(tt-1)+xp(V(tt-1))*Dt;
    %X(tt)=X(tt-1)+V(tt-1)*Dt; %same as above line
    V(tt)=V(tt-1)+vp(X(tt-1))*Dt;
end

%Try changing Dt
figure
plot(Tvec,X)
xlabel('t')
ylabel('X(t)')

ExactT=linspace(T0,T,10000);
ExactF=X0*cos(sqrt(k/m)*ExactT);
hold on
plot(ExactT,ExactF)
legend('Numerical','Exact')

%% Make a MassOnSpring Class

%see MassOnSpring.m

S1=MassOnSpring()
S1.k=1;
S1.m=1;
S1.V0=0;
S1.numInt()
S1.plotResult()

%Make another object and use different properties
S2=MassOnSpring()
S2.k=1;
S2.m=1;
S2.V0=10;
S2.numInt()
hold on
plot(S2.Tvec,S2.X,'r')
legend('V0=10','V0=1')









    
    
